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ABSTRACT 

We present a method which extends Monte Carlo studies to situations that require a large dynamic 
range in particle number. The underlying idea is that, in order to calculate the collisional evolution of 
a system, some particle interactions are more important than others and require more resolution, while 
the behavior of the less important, usually of smaller mass, particles can be considered collectively. 
In this approximation groups of identical particles, sharing the same mass and structural parameters, 
operate as one unit. The amount of grouping is determined by the zoom factor - a free parameter 
that determines on which particles the computational effort is focused. Two methods for choosing the 
zoom factors are discussed: the 'equal mass method,' in which the groups trace the mass density of 
the distribution, and the 'distribution method,' which additionally follows fluctuations in the distri- 
bution. Both methods achieve excellent correspondence with analytic solutions to the Smoluchowski 
coagulation equation. The grouping method is furthermore applied to simulations involving runaway 
kernels, where the particle interaction rate is a strong function of particle mass, and to situations 
that include catastrophic fragmentation. For the runaway simulations previous predictions for the 
decrease of the runaway timescale with the initial number of particles N are reconfirmed, extending 
N to 10 160 . Astrophysical applications include modeling of dust coagulation, planetesimal accretion, 
and the dynamical evolution of stars in large globular clusters. The proposed method is a powerful 
tool to compute the evolution of any system where the particles interact through discrete events, with 
the particle properties characterized by structural parameters. 



Subject headings: gravitation — instabilities 
protoplanetary disk 

1. INTRODUCTION 

In the field of natural sciences the dynamical evolution of 
a population is often determined by interactions that operate 
at short range, e.g., when the members come into contact, but 
which do otherwise not influence each other. The dynamical 
evolution of the system is then regulated by discrete events 
and the system can be described statistically. The most di- 
rect interaction is, of course, the merging of two bodies that 
come into contact. Examples are numerous. In the Earth's 
atmosphere water vapour condenses on cloud droplets (tiny 
water drops), which merge to form rain drops. In biology, 
the size distribution of many species of animals can be found 
by assuming animal groups merge when they make contact. 
The distribution is then determined by a power-law, which 
seems to agree with observatio ns of, e.g., schools of t una fish 
or herds of African buffaloes (Bonabe au et al.lfl999h . Sim- 
ilarly, a power-law emerges wh en bigger or ganisms are as- 
sumed to feed on smaller ones dCamacho & Sole 20011). In 
astrophysics, too, coagulation processes appear. The tran- 
sition of the dust component in protoplanetary disks - from 
(sub)micron grains to planets - covers a huge dynamic range: 
many tens of orders of magnitude in mass. Due to small 
molec ular forces, micron-sized grains will stick when they 
meet (Blum 2004). How sticking continues as particles grow 
to macroscopic proportions (~cm or larger) is still unclear but 
once the km-si ze is reached gravity will dominate the accre- 
tion process (Lissaueii ll993t iGoldreich et al. 2004). Another 
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example where coagulation is of importance are dense stel- 
lar systems. Although stellar densities are usually too low 
for collisions to be likely on a Hubble timescale, it has been 
shown that mechanisms, e.g., mass-segregation or equiparti- 
tion, c onspire to make stellar collisions feasible (Freita get all 
2006). 

The coagulation process, excluding fragmentation and in- 
jection, i s described mathema tically by the Smoluchowski 
equation (Smoluchowski 1916), which in the continuous limit 
reads 



df(m) _ 
dt 

-— I dm' 
2 J 



-f(m) 



I dm' f(m' 



)K(m, to') 



f(m — m')f(m')K(m — to', m), 



(1) 



where the terms on the RHS denote, respectively, the loss and 
gain terms of particles of mass m. The distribution function, 
f(m), gives the number density of particles (e.g., dust grains, 
stars, etc.) of mass m; i.e., 



/(to) • dTO = number density of particles in 
mass-interval (to, to + dm). 



(2) 
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Given that the initial conditions of the distribution, i.e., 
/(to, t — 0) are specified, the coagulation is fully determined 
by the coagulation kernel, K(m, to'), which is the product of 
the cross section of collision and the relative velocity of the 
particles, K(m, to') = cr(m, m')Av(m, to'). It gives the rate co- 
efficient (units: cm 3 s _1 ) for the collision between two par- 
ticles of mass to and to'. For only three distinct cases of 
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Note. — Moment expressions for the classical analytic 
kernels: Mq describes the particle density, M\ the mass 
density (which is conserved) and M% is related to the mass 
peak, m p , of the distribution. The initial distribution at 
t = is monodisperse, f(m, t = 0) = 6(m - 1). The numer- 
ical factor of 1 /2 in the sum kernel causes its evolution to 
coincide with the other kernels at small t. The expressions 
for the product kernel, K = mimj become invalid after the 
formation of the gel, at t = 1 . The mean mass (m) and peak 
mass m p are related to the moments as (m) = M\ /Mq and 
m p = Mi I M\ . 

K does equation ([T]) have a, closed form, analytic solution 
dTanaka & Nakazawall 19941) . These are: 

• The constant kernel, K — cnst.; 

• The additative (or sum) kernel, K oc (m, + mj); 

• The product kernel, K oc m^mf, 

These three kernels each show their own characteristic behav- 
ior. This can be seen by considering the moments of the dis- 
tribution 



M p (t) 






f(m, t) Am. 



(3) 



Multiplying the Smoluchowski equatio n by mP and in tegrat- 
ing over m yields the moment equation dLevv raz 2003) 



AM, 



At 



i r°° 

- = — I Am Am' K(m, m')f(m, t)f(m', t) 
2 Jo 

x [(m + m'Y - m p - m' p ] , 



p ] , (4) 

from which M p (t) can be found. In Table[T]the zeroth, first 
and second moments as function of time for monodisperse 
initial conditions are given for the three classes of kernels. 
The first moment M\, the mass density, is conserved. The 
zeroth moment, Mo, describes the (decreasing) particle den- 
sity; it therefore defines the mean mass of the population, 
(m) - M\/Mq. The second moment, Mi, is of considerable 
importance in this paper as it traces the mass distribution on 
a logarithmic scale, i.e., mf(m)Am - m 2 f(m)Alogm. We will 
therefore define the peak mass as m p = M2(t)/Mi(f), For most 
continuous distributions m p corresponds to the particle sizes 
that (together) dominate the total mass of the population. 

These three kernels are each representative of a very distinct 
growth behavior. In the constant kernel growth is orderly: 
not only do (m) and m p evolve linearly with time, their ratio, 
which indicates the (relative) spread in the distribution, stays 
constant as well. This relative narrowness of the distribution 
means that this kernel can be relatively easily simulated by nu- 
merical methods. The situation for the sum kernel, however, 
is different. Initially, (t <k 1) it evolves similarly as the con- 
stant kernel, but after f — 1 it starts to grow exponentially. The 
exponential nature of the growth not only holds for (m) or m p 
but also for the spread in the distribution, m p /(m) ~ exp[|f] 
(see Fig. [T}. 

The third class of kernels are the runaway kernels, of which 
K/j = m,mj is its analytic representative. From Table Q] a pe- 
culiarity can be noted: the moment expressions do not make 



sense after a certain amount of time: the mean density will 
become negative after t = 2, whereas Mt(t) becomes infinite 
at t — 1 . The reason for this behavior is that at t = la run- 
away particle of mass m«(?) is created that separates from the 
rest of the distribution. In theoretical studies of coagulation 
this particle is called the 'gel' and /««(£) becomes infinite at 
precisely t - 1 . In physical situations, where the total mass 
in the system is finite, m«(?) grows continuously with time 
but nevertheless separates from the distribution. This discon- 
tinuity in the distribution is also the reason why the moment 
expressions become invalid after t — 1 . 

In systems where small number statistics (one runaway 
particle) become important the assumptions underlying the 
Smoluchowski equation become invalid. Coagulation is more 
accurately described by the stochastic coagulation equation 
which assigns a probability that the system is in a specific 
state {«£}(?), with {n^} the discrete particle numbers for the 
mass bins k. The Smoluchowski equation only follows when 
certain assumptions are made, like (n,nj) = (»,)(» ,), which 
break down when i and j are statistically correlated (Gillespie 
1975). Clearly equation (Q]) is dangerous when it is used to 
describe the coagulation process - a process that is inherently 
discrete. 

In many situations of physical interest K = crAu cannot be 
expressed in the simple forms above where analytic solutions 
are available. If K can be expressed in terms of a simple mass- 
scaling behavior of the form K oc (mass/, then semi-analytic 
solutions are available (Silk & Takahashl l 19791) . However, in 
many physical situations [5 will not be constant over the whole 
mass range under consideration. For example, in studies of 
dust coagulation in protoplanetary disks, relative velocities 
can be divided into sev eral 'regimes,' each w ith a different 
velocity behavior (e.g.. lOrmel & Cuzzill2007h . Another ex- 
ample, is the enhancement of the collisional cross section by 
gravitational focussing, i.e., 

2G(m,- + mj) 



TlRt 



1 + 



(5) 



R s (Av u Y 

where R s is the sum of the radii of the bodies and G Newton's 
gravitational constant. Equation (0 displays a natural break 
in the scaling of cr with mass: i.e., <x oc m 1 ^ at low masses 
vs. cr oc ra 4 / 3 at high masses. The break happens at the point 
where the relative velocity of the particles starts to become 
less than the escape velocity u esc = -y/2G(w,- + mj)/R s of the 
biggest body, Auy/Desc < 1. 

So in most situations a numerical approach to the coagu- 
lation problem is required. Most directly, equation (fl]i can 
be integrated numerically. This entails dividing up the mass- 
axis into discrete bins and applying equation dT) to all of 
these bins. However, this is a rather elaborate procedure 
as, e.g., mass conservation is not intrinsically conserved, and 
most works have actually avoided to explicitly us e the Smolu- 
chowski equation . For example, ISpaute et alj d!991l) and 
llnaba et aUd 1999b use a (fixed) binning method that assumes 
the mass inside a bin is distributed as a power-law; coag- 
ulation is then an interch ange of mass between the various 
bins. Recently, Bra uer et all d2007l) have modified a matrix 
metho d, the Podolak algorithm, (see also iKovetz & Ofundl 
1969) which distributes the mass of a c oagulation even t 
eve nly between i ts tw o nearest grid points. Wetherill ( 1990) 
and llnaba et al.1 (1999) use a method where the mass is di- 
vided over a discrete number of 'batches,' each characterized 
by a single mass. Interactions with other batches and with 
themselves cause the batches to grow, i.e., move up the mass 



High dynamic range Monte Carlo 



10" 



10 



:l I III I I III I I III I I III I I III I I III I I III I I 1+ 



m p = {m 2 )/{m) = exp[r]| / 




il i i ill i i "I i i ill i i nl i i ill i i nl 



10° 10 1 10 2 10 3 10 4 10 5 10 6 10 7 
mass, m 



10 s 



Fig. 1. — The m f(m) mass function of the sum kernel (K/j = 4(m,- + mi)) 
at t = 15 (solid curve). At this time the mass function has achieved self- 
similarity. The mean mass, (m), and peak mass, m„, are indicated. Dashed 
lines of logarithmic slopes and 1 , respectively, trace areas of equal mass- 
and number density per logarithmic unit of mass. Due to the f(m) oc m ' 
fall-off towards small masses, the small particles dominate by number, yet 
most of the mass resides in the large particles (around m = m„). 



scale. This representation of intrinsically discrete quantities is 
especially useful in the case of non-continuous distributions 
(runaway kernels). 

Bin size techniques are generally one-dimensional, i.e., the 
quantities that enter equation (fl} (cross sections, velocities) 
depend on mass only. This approximation may neglect vari- 
ables that play a pivotal role in the coagulat ion process. For 
example, a charge distribution on grains (Marshall & Cuzzi 
2001; Konopka et al. 2005) is difficult to include in this way. 
A more subtle example concerns the coagulation of dust par- 
ticles in disks. Here, collisions affect the internal structure of 
the particles such that their porosity is not constant, which af- 
fects their aerodynamic properties and hence their collisional 
evolution. Ossenkopf ( 1993) has extended the Smoluchowski 
formalism to solve the two-dimensional distribution function 
for the collisional evolution of dust aggregates, using their 
flufhness as the second variable. However, it is clear that the 
increased number of connections between the bins, caused by 
the higher dimensionality of the problem, makes the proce- 
dure more elaborate, and slower - impractical, perhaps, for 
three or more dimensions. 

In these cases one has to resort to Monte Carlo (MC) sim- 
ulations. A good description of how a direct-simulation MC 
method is implemented is outlined by Gillespie ( 1975) which 
we will briefly summarize below (§ 12. 11 1. In MC methods 
K(m, m') determines the probability an event (here: the colli- 
sion between the particles) takes place. For multi-dimensional 
models it is probably the preferred method. Besides, the MC 
method as described by Gillespie ( 1975) provides an exact de- 
scription of the stochastic coagulation process, and will also 
accurately describe runaway kernels. However, in MC simu- 
lations one is limited to the finite number of particles that can 
be followed. It is therefore best suited for distributions that re- 
main relatively narrow, such that a modest number of particles 
gives a good knowledge of the distribution function /, On the 
other hand, when the distribution is wide a large number of 
particles must be present in order to obtain a good characteri- 



zation of /; the required particle numbers then severely limit 
the applicability of the code. 

These effects can best be illustrated in the case of the sum 
kernel. In Fig. Q] w e have plotted the theor etical value of 
m 2 f(m) at t — 15 (Tanaka & Nakazawa 1994) correspond- 
ing to Kjj - \{mi + m.j), and an initial monodisperse distri- 
bution of mass ;«o - 1 an d f( m > f = 0) = 6(m - mo). At 
t — 15 the distribution has reached self-similarity and the 
shape is the characteristic one for the sum kernel. In Fig.[T] 
m 2 f(m) is plotted showing (per logarithmic mass interval) the 
mass density distribution. Due to the m~ 3 ^ 2 fall-off of f(m) 
with mass in the low-m tail of the distribution, the smaller 
(m ~ ;«o) particles dominate by number, although it is the 
larger particles that contain most of the mass. If, for example, 
the total simulated mass is M tot , there are ~ M tol /m p particles 
representing the mass peak, and ~ M to t/(m) particles in total. 
The fraction of m p particles then decreases with increasing 
m p as {m)/m p oc nip . This behavior inevitably causes com- 
putational problems since a few massive particles remain in 
comparison to a huge number of small particles, while both 
regimes must be included in order to obtain a good charac- 
terization of f(yri). With increasing m p , a point is reached 
at which there are simply not enough particles to resolve the 
complete distribution. 

In short, the dynamic range for these kinds of simulations is 
too large to tackle with conventional MC-methods; more par- 
ticles are required but this is impossible due to computational 
constraints. To overcome these problems - i.e., to extend MC- 
methods to situations that require high dynamic range - we 
will introduce the assumption that one simulation particle can 
represent multiple physical particles: the grouping method. 
§ l2.2l introduces the grouping method. First the framework of 
a MC simulation is described in § 12.11 then § l2.21 § l2~4l out- 
line the grouping method in detail. In §[3] the new method is 
tested against the analytic solution of the sum kernel (§ 13. Il l 
and the product kernel (§ I3.21 i. In § 13.31 we consider a class 
of runaway kernels and test the assumption that the runaway 
timescale goes down with the size of the simulation (keeping 
the physical quantities such as the initial mass density con- 
stant). Next, in § !3.4l we present steady-state distributions that 
result from fragmentation of the largest particles. We briefly 
summarize our findings and discuss astrophysical applications 
in §13 

2. THE GROUPING ALGORITHM 

Before introducing the grouping method (§ 12.21) we first 
briefly outline the main features of a direct simulation Monte 
Carlo code. 

2.1. The Monte Carlo code 

The direct simulation M onte Carlo meth od in its simplest 
form consists of three steps (Gilles pigll975h : 

(i) calculation/update of the collision rates, {C,y}; 

(ii) selection the collision by random numbers: which par- 
ticles collide and the timestep involved; 

(iii) execution of the collision: removal of the collision part- 
ners and addition of the new particle(s). 

A Monte Carlo simulation starts with the calculation of the 
collision rates between the particles, C, ; = KjjfV where "V is 
the simulated volume and 1 < i, j < N p . Here, N p is the num- 
ber of particles in the simulation and Cy, the collision rate, 
is related to the probability of collision, such that Cy Af is the 
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TABLE 2 
List of symbols 



symbol 



description 



Ami 
Cu 



M k 
M u ,Mli 

N e 

N„ 
N, 
N' x 
N 
'V 
f* 



(in) 



mass dependence of kernels as in K oc rtf 

mass change of i'-particle due to accretion of y'-particles 

collision rate [s _1 ] between species i and j (including 

duplicates) 

collision rate [s _1 ] between particle i and j (excluding 

duplicates) 

coagulation kernel or rate coefficient [cm 3 s ] (velocity X 

cross-section) 

k-th moment of distribution 

first, second most massive particle 

group splitting factor 

number of groups, 2-_fj w, 

number of total particles in simulations 

number of species 

target number of species (N* ) or groups (N* ) 

initial number of particles, N p (0) 

simulation volume 

maximum allowed fractional mass change of a particle during 

group collisions 

occupancy number of species ('; number of duplicates 

characteristic mass of distribution (determining {z,} in equal- 
mass method) 

mean mass of the distribution M\ /Mq 

mass of species i 
m p peak mass of the distribution Mi jM\ 
Wj group number; number of groups for species i 
Zj zoom number of species i 

Note. — Definitions for key parameters used in this paper. 

probability that particles i and j collide in the next infinitesi- 
mal time At. In total there are N P (N P - l)/2 collision rates. A 
trick which saves memory requirements (although it increases 

CPU-time) is to store only the partial rates, C; = £ , '' Cy. The 

total collision rate is C to t = £■ '' C,\ In step (ii) three ran- 
dom deviates (7) determine the time interval to the next col- 
lision, and the particles involved (i and j). As described by 
Gillespie (1975) the time increment At is found from the to- 
tal collision rate, At = -C t ~{ In 7, whereas the particles that 
are involved are found from sampling the {C,} and {C/j}. Step 
(iii) computes the new properties of the particles that result 
from the collision. The fact that in MC-methods these steps 
are separated shows its advantages when introducing multiple 
variables; by including multiple variables, the outline of the 
program, i.e., steps (i)-(iii), will remain the same. Finally, the 
(N p - 1) collision rates associated with the newly created par- 
ticle need to be re-calculated (step (i)), after which the {C,} 
are updated. The steps then repeat themselves in a new cycle. 

2.1.1. Constant- V vs constant- N simulations 

Because each collision that results in coagulation reduces 
the particle number by one, all particles will have coalesced 
into one body after N p steps. This is the constant-^V algo- 
rithm: the simulated volume *V and total mass stays the same. 
However, in many cases this is not the intended behavior: the 
volume *V occupied by the initial number of particles in the 
simulation, N - N p (0), is only a small sample of the total 
volume under consideration. Although *V may start out as a 
volume large enough to represent the particle population, af- 
ter sufficient coagulation (and decreasing particle numbers) it 
looses this property: a larger volume is needed in order to get 
the desired sampling. However, increasing the initial num- 
ber of particles N is no solution since the number of collision 
rates becomes quickly impractical. 



A natural solution to prevent the statistical deterioration due 
to the decreasing particle numbers is to keep the numbers con- 
stant: the constant-Af algorithm. Here, the number of (simu- 
lated) particles is artificially kept constant by duplicating a 
particle after each coagulation event. Thus, mass is added to 
the system, which physically translates into an expansion of 
the simulated volume, because the volume-densit y stays the 
same. As described in Smith & Matsoukas (1998) the new 
particle is taken randomly from the existi ng (simulated) distri- 
bution and is therefore called a duplicate. Smith & Matsoukas 
( 1998) have shown that the error introduced due to the dupli- 
cation scales logarithmically with (m)(f)/mo, which is a mea- 
sure of the amount of growth. This is much improved com- 
pared to the constant- ^ method in which the reduced num- 
ber of particles causes the error to scale as the square-root 
of the growth. In principle, this method allows for indefinite 
growth. Fundamentally, the assumption is that the simulated 
volume (containing N p particles) at all times stays represen- 
tative of the parent distri bution (containing » N„ particles). 
We will test whether the Smith & Matsoukas ( 1998) duplica- 
tion mechanism also holds in systems where the growth oc- 
curs over many orders of magnitude. 

On the other hand, in situations where runaway growth oc- 
curs, the collisional evolution will depend on the initial num- 
ber of particles present in the system, even if locally the same 
physic al conditions apply, i.e., for th e same density (Wetherill 
1990; Malyshkin & Goodmanll2001l) . Clearly, a constant-AT 
simulation is not the proper way to treat these kernels; we will 
always use a constant-'V algorithm when treating these ker- 
nels. It is challenging to test simulations in which N - N p {0) 
becomes truly astronomical. Here, too, the grouping method 
will be applied. 

2.1.2. The species formalism 

A natural way to deal with the existence of a large number 
of identical particles, e.g., resulting from monodisperse initial 
conditions or through the duplication mechanism in constant- 
N simulations, is to combine them in the calculation of the 
collision rates (Soouge 1985; Malvshkin & Goodman 2001; 
Laure nzi et alj|2002c lLaurenzi & D iamond 2003). Since the 
collision rates of the duplicates are the same, it is efficient to 
introduce a multiplicity array, g, for all distinct particles. This 
is like saying we have 1 particle of type 1, 5 particles of type 2 
. . . gi of species i. It causes the number of distinct simulation 
particles (N s , or number of species) to become less than the 
total number of physical particles in the simulation (N p ) since 



N p = 



z*- 



(6) 



This algorithm is much more efficient, for it is the number 
of distinct particles (N s ) that matters in the calculation of the 
collision rates. The collision rate for one collision between 
species i and j then becomes 

Aij = gigfij (i * j); A u = \gi(gi - \)Cu (i = j), 

(7) 
where the Cy are the individual collision rates between two 
single particles. So the Aij give the probability that a (single) 
particle of species i collides with one of species j. After the 
collision partners are determined their multiplicity is first de- 
creased (g t — > gi — 1; g j — > gj - 1). If this means that g, 
becomes zero, no particles of type i exist anymore and it is 
consequently removed (A^ — > N s - 1). The new particles that 
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result out of the collision are added as a new (and distinct) 
species (N s —> N s + 1). When particles are characterized by 
only one parameter (mass) it would in this case be more ef- 
ficient to assess first whether the species is really unique, but 
in this study we treat a general approach. Duplication causes 
the occupancy numb er of the duplicated sp ecies k to increase 
by 1, gk — > gk + 1- lLaurenzi et alj d2002l) describe in-depth 
how the species formalism can be implemented and quantify 
its computational efficiency. 

We note that this grouping of the collision rates is nothing 
more than a computational trick, which, within the assump- 
tions of, e.g., the duplication mechanism does not approxi- 
mate the MC-method. Occupation numbers may become very 
large, though. For example, for monodisperse initial condi- 
tions g\ - N p initially but g\ may increase if the constant-iV 
algorithm is enforced and *V grows, as in the case of the sum 
kernel (Fig.Q]). In fact, the existence of a large number of 
identical particles is the basis of the (physical) grouping algo- 
rithm. 

2.1.3. Limiting N s . 

The accuracy and efficiency of Monte Carlo simulations is 
then determined by a single parameter, N s . Ideally, one wants 
to constrain N s . In the constant-^ simulations this means that 
N s should fluctuate around a (fixed) target value N*. (There 
is no need to enforce an exact 1-1 correspondence.) Here we 
list a few mechanisms that can be applied to balance N s . Op- 
erations that increase N s : 

• Duplication (g^ — > g^ + 1). Although duplication does 
not directly increase the number of species, the in- 
creased multiplicity of a species will decrease the like- 
lihood that it becomes empty of particles (<% = 0) and 
that it is removed; 

• Coagulation/Fragmentation. Because we generally as- 
sume that the collision products are different species, 
the number of species always increases with coaguala- 
tion and fragmentation. Especially fragmentation may 
lead to a proliferation of the number of species. 

And operations that decrease the N s parameter: 

• (Physical) grouping (discussed below). This has the 
effect that many particles are involved in the collision 
and therefore increases the likelihood of a disappearing 
species (gr,- — » 0). 

• Coagulation. Requires two particles, with which again 
a species may disappear. 

• Particle removal. This is the reverse of duplication. If 
the number of species keeps increasing, we may ran- 
domly remove particles to enhance the probability of 
obtaining <?,- = 0. 

• Merging. This is the most drastic approach, in which 
different species are merged together in a single species 
whose properties reflect the properties of its progeni- 
tors. Because of the averaging, this procedure is some- 
what against the spirit of the MC-methods. It is in fact a 
smoothening where particles of different properties are 
averaged, including their structural parameters. How- 
ever, the mean parameters may still be applicable to de- 
scribe the present state. 

In this work the species-averaging is performed as follows. If 
a species needs to be merged, we first look for its neighbor 



species in terms of mass, and then mass-average over their 
(structural) properties. As in this work mass is the only prop- 
erty, it simply means that for species i and j their merged 
occupation number becomes g — gi + g y and the new mass 
(gitrii + gjinj)/(gi + gj). (If needed, the latter quantity is 
rounded to integer values). If structural parameters would be 
present, however, this procedure can be simply extended to 
any structural property iff, i.e., if/ — > (g$i + g$ j)l(gi + gj), 
where if/ may, e.g., represent charges or filling factor. 

2.2. The grouping algorithm 

To simulate a coagulation process by Monte Carlo methods 
in a way to obtain growth over many orders of magnitude, two 
criteria have to be met: 

I The simulation must involve a large number of physical 
particles (Np), such that the population is a proper sam- 
ple of the continuous distribution. Physically, it means 
the spatial volume that is sampled will be large enough 
to account for all various kinds of particles. 

II The number of species (N s ) cannot become too large: 
since there are ~ N^ collision rates to be taken care of, 
the computational speed will slow down rapidly with 
increasing N s . 

To meet these conditions - high N p and low N s - the occu- 
pation numbers g, must be very large. This can, in principle, 
be achieved by duplicating particles many times. However, 
a new problem then arises. Since the total collision rate, C to t, 
scales with the total number of physical particles squared, i.e., 
C to t ~ ^Vp, the time step per collision is proportional to NZ 1 . 
The simulations then seem to freeze as the required CPU-time 
per unit simulation time blows up. In the MC-simulation the 
focus lies on the particles that dominate by number, which 
may sometimes be irrelevant from a physical point of view. 
This occurs, for example, with the sum kernel. Here, the most 
dominant collisional process - a collision between a small 
particle and a very massive one - is rather meaningless: a 
single collision neither affects the mass of the massive parti- 
cle nor the number of small particles. The ubiquity of many 
of these subleading collisions then overwhelms the more im- 
portant ones involving more massive particles, a process that 
together is very CPU intensive. This leads to the conclusion 
that, given (I) and (II), a third condition must apply 

III A simulation involving a large number of particles (N p ) 
yet with only 2 physical particles colliding per MC cy- 
cle can become very inefficient. Therefore, (small) par- 
ticles should be grouped together - i.e., considered as a 
single unit - during the collision. 

Due to the duplication mechanism there is a natural way to 
group particles together. In equation (O it was already seen 
that identical particles can be effectively absorbed in the cal- 
culation of the collision rates. We will now take this idea fur- 
ther and write the gr,'s (i.e., the occupation number of species 
as 

gi = Wj2 Zi , (8) 

where wi is the group number and z, the zoom number (an 
integer equal to zero or larger) of species i. Thus, instead of 
tracking g-, particles, we now simulate only w, groups, each 
containing 2 Zi particles. It is then the number of groups (N g ) 
that determines the collision rate per cycle. Because collisions 
are now between two groups of particles, involving 2 Zi and 
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Fig. 2. — Illustration of collisions in the grouping algorithm. (A, top) Two 
groups with z; = 3 and Zi = 1, respectively, collide. The zoom factors (z,) 
give the number of particles within the group (2 Z >). Less massive particles 
have higher zoom factors than particles of higher mass. All particles within 
the group are identical, (bottom) Assuming sticking, the collision product is 
obtained by equal distribution of the collision partners. The end-product is 
a new species consisting of one group where the zoom factor corresponds to 
the lowest z of the progenitor particles. (B) Same, but for in-group collisions. 

2 Z ' physical particles, coagulation is speeded up significantly. 
That is, 

N, N, 

(9) 



N o = 2 m Np = X Wi2Z '' 



and we can, with a limited number of groups (A^), repre- 
sent large numbers of physical particles (N p ) by choosing ap- 
propriately high zoom factors z,-. This satisfies (I) and (II) 
and since each collision is now between two groups involv- 
ing many particles we also avoid the inefficiency described in 
(III). Note that the base factor of 2 in equation dS) is arbitrary; 
however, it must be an integer (gr, is integer) and 2 is probably 
most convenient. More critical are the choices for the zoom 
factors, Zi'- these must be chosen to reflect the nature of the 
system that is under consideration (see § !2.3l l. 

A group collision is illustrated in Fig.|2^. The first group is 
magnified three times (zj - 3), and the second has z; = 1 
(zj > Zi in this and all other examples). The groups thus 
consist, respectively, of 8 and 2 identical particles. During 
the collision, assuming sticking, the particles are equally dis- 
tributed over the groups so that we end up with one group 
consisting of two identical particles. The mass of the new 
particle then becomes ;«,■ +4my and the zoom factor z - Zi - 1 
(the lowest of the two; the grouping may always be adjusted if 
needed). Note that in Fig. [2^ only particles of different groups 
collide: out-of group collisions. However, if their An do not 
vanish, particles within a group can also collide: in-group col- 
lisions. This procedure is illustrated in Fig. [2b. The zoom 
number decreases by one in the final configuration (group col- 
lapse). 

To generalize, for out-of-group collisions with Zj > Zi every 
/-particle collides with 2 z '~ Zi j-particles. In total there are 2 Z > 
collisions, whereas for in-group collisions (Fig.|2]3) there are 
2 Zi /2 collisions. Therefore, the collision rates of equation (0 
- the rates corresponding to one collision - are divided by 2 ZJ 
and 2 Zt ~ 1 , respectively, i.e., 



Afj = gigjCij/2 Zj = WiWj2 Zi Cij 



(Zi < zj, i * j); 
1)C«, (i = ;)(10) 



give the collision rates between two groups. Thus, although 
collisions between groups are less likely (the rate is decreased 
by a factor 2 Z >), this is 'compensated' by the 2 Z > particle- 
collisions that occur inside the group. 

However, grouping does introduce various sources of sys- 
tematic errors: 



1. the outcome of a collision between two groups is 
fixed; in reality, a distribution of configurations will be 
present; 

2. it is assumed that the grouped collision rates obey Pois- 
son statistics like the single particle collision rate; 

3. the individual collision rates {C,y } do not change during 
the group collision process to ensure that equation ( fTOb 
is applicable. 

A fixed configuration for the outcome of the grouping pro- 
cess is required because it prevents (unwanted) proliferation 
of the diversity of species, which would increase the N s pa- 
rameter and strain CPU requirements. However, the diver- 
sity of the simulation is already specified by the N s parameter 
and we may simply check how N s affects the accuracy of the 
collisional process. This point is similar to the second con- 
cern where all collisions are forced to occur at one point in 
time. Also, note that we must use the memoryless property 
of Poisson statistics; otherwise all collision rates must be ad- 
justed at each timestep. However, it can be shown that, pro- 
vided the collision rates stay constant, the grouped rates in 
equation ( TT~Ob will conserve the mean of the particle distribu- 
tion if i + j (see Appendix lAl. The last assumption is the 
weakest: the individual rates between the particles in princi- 
ple change after each individual collision as the collision af- 
fects the properties of the particles. Therefore, expressions as 
g-,gjCij should actually be averaged over the number of col- 
lisions that make up the group process. This is of course 
a very tedious procedure with which we would go back to 
the one-by-one treatment of the particles. Rather, we pre- 
fer to introduce a new parameter, / e , that limits the amount 
of mass accretion onto a particle during the group collisions, 
thereby strengthening the approximation that the Cy stay con- 
stant during the group collision process. 

Choosing the particle labels again such that j has the high 
zoom factor (zj > zi) the fractional mass accretion on a single 
i particle is Ato,/to; = 2 z '~ Zi mj/mi. From the preceding dis- 
cussion it is this quantity that should be limited from below 
by a fixed value, f e . One way to implement this is to reduce 
the y'-particles that take part in the collision. For example, 
in Fig. [2^ this could mean that only 2 j-particles are involved, 
i.e., only a fraction of the particles in a group. Let this fraction 
be denoted by 2~ Ne with N € > integer. Then, only 2 z i~ Zi ~ N ' 
/'-particles are involved in the collision and the mass of an i 
particle increases by 2 Z - 

Anij 



-z,-N e , 



2 Z > 



u 



(11) 



TO; TO; 

Solving for the group splitting factor, N e , gives 

N e = [-log 2 (f e 2 z 'mi/2 z 'mj)], (12) 

where the square brackets denote truncation to integer num- 
bers > 0. Consequently, the collision rates, instead of 
equation ( flOb . now become 

A? = gigjC u /2 z r^ = w . w j2 z ' +N <C if , Z; < Zj , (13) 

(the collision rates for in-group collisions do not change since 
Af e = 0). Another consequence is that the group numbers are 
real numbers, i.e., since a fraction 2~ Ne of the y'-group is re- 
moved, its group occupation number decreases accordingly, 
Wj — > Wj - 2~ Ne . The {g/} cannot become fractional, though, 
since this would physically split a particle. Also, a situation 
where < to,- < 1 is not allowed, because Aw - 1 for the par- 
ticle of the lowest zoom factor; when it arises we 'demagnify' 
the group (z; — > Z; - 1) such that w, > 1 (see below, § I2.41 i. 
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2.3. Choosing the zoom factors 

We have implemented the grouping method in such a way 
that the zoom factors, are free parameters, and there is in- 
deed considerable freedom to choose (an algorithm for) the 
{Zi}. However, it is also an important choice as it will de- 
termine how efficient the Monte Carlo simulation becomes. 
Here, with "efficient" it is meant that the number of species 
N s should not (suddenly) increase. A way that observes this 
constraint fairly well is to divide the mass inside the simula- 
tions equally over the groups. Additionally, we present a more 
general method, in which the zoom factors are determined by 
the particle distribution. 

2.3.1 . Equal mass approach 

In this method, we assume there is a characteristic mass 
m* after which the zoom factors are determined. Particles of 
mass larger than m* are ungrouped (z - 0), whereas particles 
of mass m < m* have zoom factors such that the group mass 
is around m* , i.e., 



[log 2 (m*/m,)] 



(14) 



In the case of constant- V simulations and monodisperse initial 
conditions of unity mass the simulation then starts with N/m* 
groups, where N is the initial number of particles. Therefore, 
N/m* is a measure for the number of species parameter N s , 
The target value of N s , N*, and N then determine the value of 
m* . For the constant-Af, simulations, where the mass in the 
simulation varies, the mass peak m p is a natural value for the 
characteristic mass, m* = m p . It forces most of the groups to 
represent particles around the mass peak, with those of higher 
mass ungrouped. 

Thus, the equal mass approach resolves the mass density 
of the distribution, instead of the particle density as in con- 
ventional MC-methods. Once species becomes insignificant 
in terms of mass, m{l Zi < m*, occupancy number approaches 
unity (wi - 1) and the species will be "removed" from the 
simulation when it collides at the next collision, completely 
independent of the actual number of particles a group con- 
tains. 

2.3.2. Distribution approach 

A different way to choose the zoom factors is to force the 
distribution to be resolved over its entire mass range. This 
can, for example, be done by dividing the distribution into 
exponentially-spaced mass bins and to choose the zoom fac- 
tors such that there are a fixed number of groups per bin. Here 
the bins are exponentially distributed between a lower mass 
(most often the initial mass, mo) and the largest mass Mu 
present in the simulation. For example, let Nu be the total 
occupancy number of bin k, then the zoom numbers are de- 
termined by the requirement that an equal number of groups 
is present in each bin, i.e., Nk/2 Zt is constant for each bin k. 
There is no constraint on the number of species per bin, but 
as this quantity is always less than the number of groups it 
is limited as well. For a smoother progression of zoom fac- 
tor with mass we allow for a linear interpolation within each 
bin, such that the zoom factors become a continuous, piece- 
wise, and decreasing function of mass, z{m). Thus, z is now 
determined, not by the particle or mass density, but by the 
sheer occurrence of particles of a certain mass. The bins and 
zoom factors are dynamically adjusted, as the simulation pro- 
gresses, and the distribution changes. 



The distribution method has the consequence that (high- 
m) fluctuations are resolved very well, even if they are com- 
pletely insignificant by number and mass. The zoom factors 
of these particles are much reduced in comparison with the 
equal mass method. For example, if the mo = 1 particles 
are still at zoom levels of zo situations where 2 Zo mo » 2 Z/ m,- 
are frequently encountered. A collision between a massive 
and a m - 1 group then takes place at a high value of N € 
(eq. US)) and collision rates involving high-ra groups increase 
(eq. |[T3l ). Since each collision creates a new species, the 
number of species (N s ) also increase rapidly - too rapid, in 
fact - and we are forced to merge species (§ 12. 1.3b . Merging 
of species is therefore a necessary feature of the distribution 
method. 

2.4. (De)magnification 

In the grouping method the zoom factors of the groups are 
determined in accordance with one of the two methods above. 
The zoom factors will therefore change during the simula- 
tion. When Zi increases (z,- — > Zi + 1) we speak of magnifi- 
cation, whereas when zi decreases the group is demagnified 
(zi —> Zi - 1). This dynamical adjustment of the zoom fac- 
tors is the main reason why the method is capable of achiev- 
ing very high dynamic range. Demagnification furthermore 
occurs when the occupancy number of a group becomes frac- 
tional, < Wi < 1 . However, we note again that the actual 
number of particles of a species, g„ is always integer. 

There is a small price to pay in terms of overhead for these 
re-adjustments. Because the {z,} are an argument of the colli- 
sion rates (eq. |[T3l ) these rates (or, more accurately, the {/!,} 
quantities, as these are the ones that are stored in the com- 
puter's memory, ~ N s in total) have to be adjusted at each 
(de)magnification event. This may slow down the simulation 
by a factor of 2. Finally, we note that (de)magnification is, 
within the context of the grouping algorithm, an exact proce- 
dure; collision rates are not approximated and the total mass 
of the simulation stays conserved. 

2.5. Grouping Summary 

We have outlined an approach for a Monte Carlo coagu- 
lation scheme that can involve many particles - allowing the 
growth to be orders of magnitude - yet avoids the severe com- 
putational problems that would naturally arise with a single 
particle per collision treatment. In Fig. [3] we summarize the 
main steps of the algorithm, within the context of the grouping 
method. At its core lies the approximation that many identical 
particles are present, due to monodisperse initial conditions or 
the duplication mechanism, which can be efficiently grouped 
together or (in the distribution method) merged. Collisions are 
now between two groups of, possibly, many particles and this 
significantly speeds up the computation. While it is clear that 
such grouping inevitably approximates the collision process, 
we have advocated a fundamental shift in the implementation 
of MC methods, i.e., to focus on the particles in a distribu- 
tion that really matter - e.g., those around m - m p - and 
we anticipate that the benefits resulting from the new empha- 
sis will more than justify the approximations inherent to the 
grouping method. In any case, to validate the new method a 
comparison must be made with the analytic solutions of the 
classical kernels, after which the algorithm can be applied to 
more physically interesting kernels. 

3. RESULTS 
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1. Setup state vector 

species group number zoom factor Particle properties 

mass (other properties) 

1 W\ Z[ mi 

2 wi zi m2 



Note: wj can be real but not less than unity (w, > 1); z, is integer. Number of 
particles of species i, g, = 2 Zi w h is also integer. 

2. Calculate group collision rates 

Afj = WiWjQj J2 Z < +N - (if!*./); 

^u - \ w t( w i - l)C,v (/ = j, In-group collision) 



Store partial rates, A, = V Af k 
Store total rate, A tot = V i, 



Where: 



- Cif individual collision rate between two particles, Ky/"V. The 
collision kernel K u = o-,j(A«),j depends on mass and other parti- 
cle properties. 

- "V: simulation volume 

- N/. group splitting factor that regulates mass accretion of j- 
particles on /-particles (here: z, > z,) to avoid excessive instan- 
taneous mass load. It is defined as 



N € = \-log 2 .(f t ym t /Wmj)\ 



d) 



where f c is the maximum allowed fractional increase in mass of 
the larger particle /. 
- N s : total number of species 

3. Pick collision particles 

This follows the "full conditioning" method of Gillespie (1975). 

time increment: At = A~~} { ln(l/r); 
i-i i 

i - particle: ^ A k < ?A m < ^ A k ; 



7-1 



./ 



-particle: £ A? k < rAf < £ 4 



where r are (different) random deviates. 



4. Perform collision 

- (Choose/rearrange /, j such that z, > z,). 

- Determine the group-splitting factor N € and the collision multi- 
plicity 

- Each /-particle collides with N = 2 : '~ Z >- N ' ^'-particles. 

Perform the group collision: the sequential collision of N ./-particles 
with the /-particle. Determine the new (structural) properties of the 
particles created. 



5. Update collision rates, state vector 



1, Wi 



- Reduce group-number of collision species: u>, -> 
wj - N e . 

- If wi or Wj becomes 0: remove this row from the state vector 
(N s — > N s - 1) and remove entry from {Af}. 

- Add the new particle to the state vector at index N s + 1 (N t -» 
N s + 1). 

- Update the {Af}. For every entry k look whether the rates A ik , A jk 
should be subtracted (if k < /, resp., k < j) and add rate A qk if 
k > q. Recalculate A m . 

6. Adjust I:,!, merge species, choose duplicates 

These steps are necessary for a smooth progression of the program, 
but do not have to be processed at every cycle. Adjust collision rates, 
{A,} accordingly. Zoom factors can be adjusted up or down: 

- magnification: zoom factor of species is increased, z, -> z, + 1 
and Wi — > Wj/2. 

- demagnification: zoom factor of species is decreases, z, -> z, - 1 
and w t — > 2wj. 

Note: We have outlined two methods to regulate fc}. (I) Equal mass method: 
z.i is related to the mass peak and requires duplication; (II) distribution method: 
z.i determined by the distribution. In both methods the number of species pa- 
rameter iV s must be constrained. 

- Merging: properties of similar species are averaged to one species. 
This limits the number of species in a mass interval. 

- Duplication (constant-/V s simulations only): pick species-^ ran- 
domly from statevector (with weights {«;,)) and add to it (w k -> 
w k + 1). Update collision rates {A,}. 

Note: Duplication does not directly increase N s but increasing group numbers 
make occurrence of empty rows less likely. If particles are duplicated the 
simulated volume "V also increases. Constant-'V simulations do not require 
duplication. 

End of cycle -> return to 3 



Fig. 3. — Chart outline the steps required to implement the grouping algorithm. 



In this section the grouping method - both the equal mass 
approach, as well as the distribution method - is tested in sit- 
uations that require a high dynamic range. Its performance is 
then discussed. First, the grouping method is tested against 
the sum kernel. We show that both the equal-mass and the 
distribution methods are capable of following the analytic dis- 
tribution, whereas conventional MC methods (without group- 
ing) will fail. The grouping method is also tested against the 
product kernel, which is one of a class of kernels that pro- 
duces runaway-growth behavior. Here, we will see that it is 
important to accurately follow the fluctuations of the mass dis- 
tribution, which can only be done by the distribution method. 
More runaway-models are discussed in § 13.31 Finally, we also 



perform simulations including a simplified form of fragmen- 
tation, where the particles are catastrophically fragmented at 
the largest scale mf rag to be injected back at the smallest scale 
niQ (§ |3.4t . We test whether the grouping method can also 
be applied in these cases. An overview of the simulations 
is given in Table [3] including the adopted numerical param- 
eters. These are chosen such that the CPU-time per sim- 
ulation is typically several hours on a modern desktop ma- 
chine (1.8 GHz clock speed). Similar CPU-times apply for 
the choice of the numerical parameters of the simulations in- 
cluding grouping. 

For the product and other runaway kernels, the system will 
after a certain time be dominated by a single particle. Hence, 
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TABLE 3 
List of simulations 



Kernel/description 


Figure ref. 


Type 


Grouping 


n;/n; 


N 


fe 


N ■ 

1 * sim 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


sum kernel 


Figure [4] 


constant-iVp 


none 


10 s 




io- 4 


40 


sum kernel 


Figure [5^ 


constant-JV, 


equal mass 


2x 10 4 




io- 4 


40 


sum kernel 


Figure [5p 


constant-'V 


distribution 


10 4 


10 ioo 


IO" 4 


40 


product kernel 


Figure [6] 


constant-'V 


equal mass 


2x 10 4 


10 20 


io- 4 


40 


product kernel 


Figure [JJ 


constant-'V 


distribution 


10 3 


10 20 


io- 4 


40 


runaway kernels 


Figure [8] 


constant-'V 


distribution 


~ 10 4 - IO 5 


30CM0 160 


IO" 4 


10-100 


sum: fragmentation 


Figure [5] 


constant-iV, 


equal mass 


2 X 10 4 




i<r J 


l a 


constant: fragmentation 


FigurellOl 


constant-iV, 


equal mass 


10 4 




io- J 


l 11 




Note. — Summary of all simulations. (1) Collision kernel. (2) Figure reference. (3) Simulation type: fixed or 
adjustable volume. (4) Grouping method: outlines in which way the zoom factors, {z,|, are chosen (no grouping, 
equal mass, or distribution method). (5) Target number of species (equal mass method), or groups (distribution 
method) in simulation. When multiple values apply the largest is given. (6) Initial number of monomers, or total 
mass in simulation for constant-'V simulations. (7) Fractional mass accretion. (8) Number of simulations performed. 
The initial distribution is monodisperse in all simulation with mo = 1 . 

a In fragmentation simulations a steady state is reached in which the distribution is averaged ~ 100 times during the 
same simulation. 

particles, N p , is kept constant at 10 5 . The distribution function 
is plotted at various times during its evolution, t - 1, 10, 20 
and 30. Times are dimensionless, i.e., in units of (,K"oo«o) 
where Kqq and «o are the kernel and number density corre- 
sponding to the initially monodisperse conditions of particles 
of mass mo — 1. In Fig.|U as well as in all other figures 
that show distributions from the MC-simulation, f(m) is de- 
termined using exponentially-spaced bins every factor of 2 in 
mass. Each point shows the average and spread of in total 40 
model runs. The analytic distribution functions at these times 
are given by the solid curves. 

Initially, r < 10, the agreement between the MC simula- 
tion and the analytic curve is very good. The mass range in 
the simulation is still rather modest and the 10 5 particles suf- 
fice to represent all relevant masses. However, while the peak 
mass progresses steadily to higher to, the majority of particles 
stays behind at low masses; that is, although both (to) and m p 
evolve exponentially with time (Fig.[T), m p /(m), which is a 
measure of the dynamic range, also increases exponentially. 
Since MC-methods sample the number density, increasingly 
fewer particles become associated with the mass peak, and the 
simulation looses numerical resolution. The effects of these 
low-number statistics are a bit alleviated by averaging over 
many different simulation runs but it is not merely a statistical 
deterioration: the mass peak dominates the collisional sweep- 
up of the system and the fluctuations caused by an inaccurate 
characterization blow up. There is especially a trend to lag the 
analytic solutions because, it is impossible to 'recover' from a 
lagging distribution, because the duplication mechanism does 
not duplicate the particles that are not present. At t — 20 these 
problems already become apparent when the spread is very 
high, especially around the mass peak, and there is a slight 
tendency to lag the analytic curve. At t — 30 the effects are 
disastrous: 10 5 particles are simply too few to sample a large 
enough volume to represent the large particles. 

Due to the duplication mechanism many duplicates are cre- 
ated and N s slowly decreases. This somewhat speeds up 
the computation at larger times since collisional rates scale 
with N 2 S . A mor e accurate result w ill be obtained by keep- 
ing N s constant dOrmel et alj 120071) . However, now N p in- 
creases and the simulation suffers computational freeze-down 
as it progresses. Fundamentally, either approach - the con- 
stant N p or constant N s - is the same since each cycle only 
treats one collision and each collision is equally important. 
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Fig. 4. — Results for the sum kernel, Kjj = \ (m, + mi), without grouping for 
N p = IO 5 . Distributions are plotted at dimensionless times t = 1, 10, 20 and 
30 and error bars show the spread over 40 simulations. Particles are binned 
at each factor of two in mass. The correspondence with the analytic results 
(solid lines) deteriorates over time. After t =* 20 the mass peak is no longer 
resolved by the simulation. 

no self-similar solutions exist and in order to study the evolu- 
tion properly, the simulation volume of these systems must be 
well defined. These simulations are therefore modeled using 
constant-volume simulations. In most cases the distribution 
method is used, which also requires a fixed simulation vol- 
ume. In the other cases, we may assume that a small volume 
is representative of the parent distributions, which can then be 
modeled using constant-number simulations. The only excep- 
tion is Fig.[5}5, where we have used the distribution method in 
order to compare it to the results of the equal-mass method 
simulations of Fig. [5^. 

3.1. Test case I: the sum kernel 

Before presenting the results of the grouping mechanism, 
a case without the grouping method is considered first, illus- 
trating some of the 'deficiencies' mentioned in §Q] We will 
not consider the constant kernel, because in the constant ker- 
nel the dynamic range (m p /(m)) stays confined and a group- 
ing method is not required. In Fig. |4]the number of simulation 
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Fig. 5. — Results for the sum kernel using the grouping algorithm. Particle distributions (symbols) are shown at times t = 20,40,60 and 80. Particles are 
binned at each factor of two in mass and results from 40 separate simulations are combined. The spread in the simulations is denoted by the error bars. (A) Equal 
mass method, N s = 20 000. Thanks to the grouping algorithm most of the groups are centered on the mass peak, which favours an accurate computation of the 
collisional evolution, without introducing systematic errors. (B) distribution method with N* = 10000, distributed over 20 bins. The zoom factors are adjusted 
to resolve the distribution over its entire range. 

Clearly, this 'socialist' approach is ineffective as Fig.|4]illus- 
trates: the more massive particles deserve more attention than 
their smaller brethren. 

Results with the grouping method applied are shown in 
Fig-El In Fig- Eh we use d the equal mass method for the 
groups (each group has a mass of ~ m* - m p if m < m p and 
particles with m > m p are not grouped) and choose the zoom 
factors accordingly (see eq. Q3)). The number of species, N s , 
is stabilized at a target number of N* - 2 x 10 4 through the 
duplication mechanism. That is, if N s due to coagulation falls 
below N*, duplicates are added to the system until N s > N*. 
The simulation total mass and volume then increase, but their 
ratio (the mass density) is fixed. As m p becomes larger during 
the simulation, zoom factors of particles that do not take part 
in collisions will increase, according to equation ( fl4l i. This 
'magnification' causes the occupancy numbers w, to stay low, 
despite the ongoing duplication events. When Wi reaches 1, 
the species will disappear at the next collision, freeing up this 
slot for more massive particles. The grouping thus forces the 
particles to act collectively. This procedure prevents the sys- 
tematic errors (N p constant, Fig.|4| or computational freeze 
down (constant N s ) of the ' 1 collision per cycle' approach of 
the conventional MC-methods. 

FigureEh shows that the grouping mechanism resolves the 
mass peak well during its entire simulation. The agreement 
between the analytic and numerical method remains very 
good during the entire simulation. FigureEh proves that once 
the mass peak is resolved, the simulation will accurately com- 
pute the collisional evolution, and systematic errors as in 
Fig. [4] are prevented. Thus, the resolution of the mass peak 
is key. However, at lower mass densities the groups disappear 
from the simulation as they do not have sufficient mass. 

These low-m statistics can be accurately followed when we 
switch the method for assigning the group's zoom factor {z,} 
to the distribution method, in which the zoom factors are ad- 
justed such that there are an equal number of groups per ex- 
ponential mass interval. The zoom number for the low-m par- 
ticles is then much lower in comparison with the equal mass 



method. Figure|5p shows that the mass density as well as 
number density are accurately followed over the entire distri- 
bution. It can be seen in Fig.Eb that the match to the high-m 
tail slightly deteriorates over time. This is a consequence of 
the equal spacing in log space of the bins that determine the 
zoom factors. The high-m tail, which falls off exponentially, 
becomes associated with only a single bin, and is then less 
well resolved. A more flexible placement of the bins is ex- 
pected to give an even better correspondance. Note again that 
the bins are only needed to determine the zoom factors; there's 
no actual binning of the particles in the MC-simulation. 

Although Fig.Ep more precisely follows the distribution - 
i.e., it both follows the mass density as well as the number 
density - it requires to merge species in order to restrict their 
numbers (see § 12.1.31 ). For these test cases, in which only 
mass is a parameter, Fig.Ep shows that this is not problem- 
atic. However, when structural parameters are involved, the 
merging procedure is less obvious; i.e., are species of com- 
parable mass merged or, say, species of comparable poros- 
ity? The merging principle is, therefore, somewhat against the 
spirit of MC-simulations, as there is no physical equivalent of 
'merging' in nature. Thus, for a MC simulation with many 
structural parameters, where it is not required to resolve the 
number density, Fig. Eh is probably the preferred method. In 
the other test cases of this section we will therefore first pur- 
sue whether these can be modeled by the equal-mass method. 

3.2. Test case II: the product kernel 

The kernel K,j - nijinj is the only runaway kernel with an 
analytic solution. In models that show runaway-behavior a 
particle will appear that separates from the continuous distri- 
bution to consume, in a relatively short time, all the available 
mass. This separation is an intrinsic feature of the kernel and 
is not caused by a poor sampling of the distribution. The point 
at which this separation occurs, f«, depends on the properties 
of the kernel and, in some cases, on the initial number of par- 
ticles (AT). In the case of the product kernel it always occurs 
at t - 1 , independent of N. The runaway growth properties 
of other kernels are discussed in § 13.31 
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Fig. 6. — Results for the product kernel, Kjj = mjinj, with N = 10 . The 
simulated distribution is plotted at times t = 0.4, 0.7, 0.95, 1.0 (from left to 
right, crosses) and t = 1.05, 1,3 and 2.0 (from right to left, open circles). 
The analytic solution is plotted by solid curves for t < 1 and by the dashed 
curves for t > 1 (the t = 0.95 and t = 1.05 curves almost overlap). Error bars 
denote the spread in the 40 simulations. After (» la runaway particle of 
mass m ~ N appears that consumes most of the mass in the system. 

In § [1] it was noted that the Smoluchowski coagulation 
equation becomes invalid in the case of the runaway kernel: 
it keeps assuming that the distribution is continuous and it 
is therefore ill-suited to describe a discrete system caused by 
the runaway particle. Wetherill ( 1990) has adjusted the origi- 
nal (incorrect) solution of Trubnnikov ( 1971) by inclusion of 
a term for the runaway body in equation |[T). Although not 
mathematically rigorous, this approach captures the physics 
of the runaway: as soon as the (incorrect) numerical solution 
exhibits a f(m) °c to~ 5 ^ 2 power-law tail, the interaction be- 
tween the higher mass bins causes the runaway growth. In 
the Monte Carlo simulations we also expect that the appear- 
ance of a power-law tail is a signature that runaway growth is 
incipient. 

In Fig. [6] the distributions corresponding to the product 
kernel with N - 10 20 are plotted at dimensionless times 
t = 0.4,0.7,0.95, 1.0, 1.05, 1.3 and 2.0, using the equal mass 
method. Since we start the simulation already with a huge 
number of particles, large zoom numbers are required. We 
have fixed the mass peak at to* = A/'/lO 4 , which causes 
the particles to be sufficiently magnified, except for the run- 
away particle. (The 10 4 factor may be thought of as a fudge 
parameter akin to N s or N p in the c onstant- N sim ulations.) 
Theoretical curves correspon d to the [Trubnnikovl (119711) so- 
lutions which, as Wetherill (1990) showed, are still valid for 
the continuous part of the distribution. At t - 1, runaway 
growth is expected and the continuous power-law distribution 
reaches its furthest extent: the -5/2 power-law exponent. Af- 
ter the runaway particle has separated the continuous distri- 
bution shrinks, especially in the —x direction since the more 
massive particles are consumed rapidly by the runaway body, 
which itself is not shown in Fig. [6] 

The runaway growth is generally well modeled by the 
Monte Carlo method. The distribution becomes most sen- 
sitive near t — 1 where a small change in time means a 
big change in the high-TO distribution function. Due to the 
stochastic behavior of the Monte Carlo method, the precise 
timing of the runaway event does not exactly coincide with 
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Fig. 7. — Statistics for the simulations involving the product kernel. As 
function of the mass of the largest species, Ml\ , are shown: (i) the mass of 
the second-largest species, M/,2, according to the equal-mass method (thick 
solid curve) and the distribution method (thin solid curve); and, (ii) the time 
(dashed cur\>e). Values are averaged over 40 runs. The auxiliary lines Ml2 = 
Mu and M12 = TV 2 ' 3 are also given (dotted lines). 

t = 1. Therefore, the t — 1 distribution lags the theoretical 
curve since at t — 1 either runaway growth did not yet occur 
(and the distribution still has to reach the power-law curve) or 
has already happened (and the distribution is on its way back 
from the power-law curve). In Fig.UJthe relation between the 
two most massive groups is denoted by the thick solid line, to- 
gether with the simulation time t as function of Mu . Figure [JJ 
shows that until t — 1 there is little difference between the 
mass of the most massive particle - actually group - Mu and 
the second-most massive group Mui. However, as runaway is 
reached at Mu ~ 10 5 the masses diverge very quickly. 

These results, however, are in disagreeme nt with the theo- 
retical study of Tanaka & Nakazawa ( 1994). Comparing the 
solutions to the stochastic and the statistical equation they find 
that, in the case of the product kernel, the statistical equa- 
tion (eq. HI) becomes invalid after Mu ~ N 2 ^ . It is only at 
this point that Mu and M12 should start to diverge. How- 
ever, in the equal-mass method (Fig. [71 thick line) the di- 
vergence happens much earlier: M L i never reaches beyond 
10 8 «: (10 20 ) 2 ^ 3 . This discrepancy can be explained as an ar- 
tifact of the equal mass method. Because each group is still 
of considerable mass, to* ~ 10 16 , runaway growth is just the 
collapse of one group by a continuous series of in-group col- 
lisions (Fig.UJ}); that is, until z, = 0, the Mu group does not 
correspond to a single particle and therefore never really sep- 
arates. It is thus not a proper modeling of a runaway process. 
In order to accurately follow the behavior of the particles dur- 
ing runaway growth, the groups should be demagnified much 
more rapidly than equation ( fT4l ) allows. It is perhaps sur- 
prising that the equal-mass method still works fine - i.e., its 
distribution agrees with the analytical case (Fig. [6]) - despite 
the fact that it does not resolve the runaway particles. How- 
ever, in general this is not the case: the fluctuations present at 
high-TO require more resolution than what is attributed to them 
through the equal representation of mass density (eq. |[T4l ). 

In Fig. UJ the Mu,Mi2 statistics of the product kernel are 
recalculated with the distribution method for the zoom factors 
(thin solid line). According to the Tanaka & Nakazawa ( 1994) 
prediction Mu and M12 should then separate at Mu ~ N 2 ^. 
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TABLE 4 
Runaway models 



type conditions 



growth mode 



(i) v < l,fl < 1 orderly growth 

(ii) v < I,j8 > 1 runaway growth at (o 

(iii) v > 1 instantaneous runaway growth" 

Note. — Properties of runaway models according to 
the study of Lee (2000). The exponents v and fi give the 
mass dependence of the kernel on a heavy and light par- 
ticle, respectively, with /? = v + // the total mass depen- 
dence, Kjj oc m^trt. (rrtj 3> mj). 
a Instantaneous in the limit N — » oo. 

This is exactly the trend seen in Fig. [7] We therefore con- 
clude that the distribution method is the appropriate choice to 
compute the collisional evolution of runaway kernels. 

3.3. Strong Runaway kernels 

There is an extensive literature on the phenomenon of run- 
away growth, especially mathematical (N — > oo), where the 
phenomenon is termed gelation. In general, three growth 
modes can be considered (see Table 0): (i) kernels where run- 
away growth does not appear, as in the constant or sum kernel; 
(ii) kernels that exhibit growth after some characteristic time, 
as in the product kernel; and (iii) kernels in which gelation 
is instantaneous 1 . The occurrence of instantaneous gelation 
for the class of kernels K cc (minijY with v > 1 was con- 
jectured by Spouge dl985l) on t he basis of MC-coagulat ion 
experiments and was proved by I Jeonl d!998l) . iLed ([2000) il- 
lustrates these distinct modes of solution by a simple coag- 
ulation model, in which the distribution is approximated by 
a characteristic mass m, e.g., the peak mass, and a test parti- 
cle of mass M > m. Particles of mass m dominate the total 
mass and hence n(m) oc m . Furthermore, the coagulation 
kernel is characterized by two exponents, v and id, that give 
the dependence of K(nti, mj) in the case of a heavy (/) and a 
light particle (j), i.e., K,j oc m*W, for ;«,■ » nij. The total 
mass dependence is then given by (3 — fi + v. For example, 
in the case of the product kernel v = /j. — 1 and /3 - 2. Set- 
ting the physical density equal to unity, the equations for the 
simplified merger problem become (Lee 2001) 



— = w- 
d* 2 ' 

dM „ „ 

= nfM v . 

At 

These can also be combined to obtain 



AM _ IMV 
dm \m) 



1/(1-/3) 



Equation ( [15] ) can be solved as 

[1 + 3(1-0] 

exp[if] 

where we have for simplicity taken mo - m(t 
Equation ( fTTI i can be solved in terms of m as 



for/3* 1; 
for/3 = 1, 



M(m) = i ( M o " 
I M m 2 



2m l 



-2) 



l/(l-v) 



forv* 1; 
for v = 1 , 



(15) 
(16) 

(17) 

(18) 
0) = 1. 

(19) 



1 We note that in a physical situation where runaway growth tends towards 
instantaneity, causality should not be violated. As we will see, N, even in 
astrophysical circumstances, does not become so large that the process is 
instantaneous. 



where Mo > 1 is the inital mass of the test particle. Three 
qualitatively distinct growth modes follow from these equa- 
tions: (i) if v < 1 and ft < 1 both equation (fT8l and 
equation ( fl9] l have finite solutions at all t: orderly growth, 
(ii) If v < 1 and (5 > 1, equation (TT8l reaches a singularity 
at a time t — t R — 2/(J3 - 1) > 0, 2 while equation (fT9l does 
not: both m and M become infinite at the same time f« that 
does not depend on Mo- (iii) If v > 1, equation ( fT9l reaches 
infinity for a finite value of m. Furthermore, the gelation time 
decreases with increasing Mo, ?r - Ml~ v /(v - 1). Then, if 
the initial mass in the system goes to infinity, a sufficiently 
strong fluctuation (large Mo) will always be present to satisfy 
the singularity condition for M(t) even as t — » 0; this is the 
physical reason behind the instantaneous gelation conjecture. 
Here, we will concentrate on the case v > 1 and test the 
instantaneous gelation conjecture for kerne ls K - (mimj? 
with v — 2, 3. Malyshkin & Goodman (2001) extend the sim- 
plified merger model discussed above to a model with only 
two kinds of particles - one predator and many preys - and 
one collisional process (the predator feeds of preys). This 
is the monotrophic model, e.g., one particle ('mono') that 
nourishes ('trophein') o f the o thers. Considering this model 
Malyshkin & Goodman (2001) compute the probability that 
the predator occupies any finite mass, which decreases with 
increasing time. Consequently, the gelation time f« is com- 
puted as 



t R oc (log N) 



(20) 



where N is the initial number of particles present. Thus, the 
timescale for runaway growth decreases as function of the ini- 
tial number of particles in the system, but only logarithmi- 
cally. Equation ( f20b may be anticipated from the discussion 
above on the simple merger model. Consider a distribution of 
particle masses N(m) at a fixed time. It is quite likely to expect 
that the tail of the distribution is exponentially suppressed as 
in the case with Poisson statistics, i.e., N(m) oc Nexp[-m] 
with Af the total number of particles. Setting jV = 1 we 
thus obtain Mq ~ n ~ log AT. Hence, the logarithmic de- 
pendence on the number of particles. iMalyshkin & Goodman! 
(2001) test this model by Monte Carlo simulations and found 
qualitative agreement: the gelation time decreases as function 
of log N, although the power-law exponent is not precisely 
that of equation d20b: 1 — v. However, the simulations of 
Malvshkin & Goodman (2001) were limited in their dynamic 
range; we will extend their study to much larger values of N 
to see whether equation d20b holds over a sizeable range of 
log AT. 

Note that as conjectured earlier, a good census of the fluc- 
tuations is the key to an accurate modeling of the colli- 
sional evolution. Therefore, the grouping must then be im- 
plemented by the distribution method. There is, however, an- 
other catch. At t — the number of groups of the monodis- 
perse species is wq ~ N* <k Af and the collision rate for a 
collision between two monomers according to equation (TTOt 
is An ~ NwqCh ~ wq. Then, the timescale for collisions 
involving monomers is t\\ - A^ ~ wZ l , which can, however, 
become comparable to, or larger than, the runaway timescale 
of equation ( 1201 . Therefore, the monodisperse groups have to 
be demagnified such that 1 1 \ <k t R . By forcing the appropriate 
collision timescale for the monodisperse collisions, the fluc- 
tuations are adequately resolved, but it also means that the 

2 Note that this expression disagrees with the product kernel in which 8 = 
1: the simplified merger model does not give the correct prefactor, see ILed 
(2001). 
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Fig. 8. — The runaway time t^ as function of the initial number of particles 
N for the kernel Ku = {mjmj) y with v = 2 (solid line) and v = 3 (dashed line). 
Simulations are run at N = 300, 10 5 , 10 10 , 10 20 , 10 80 and 10 160 . The power- 
law index is indicated at every interval of log N. The agreement with the 
theoretically expected value for the power-law index, 1 - v, becomes better at 
higher N. At N = 10 160 a simulation at a higher numerical resolution gives 
an even steeper slope (indicated below the line) confirming the converging 
trend towards, respectively, -1 and —2. 

CPU-time increases, especially at large N. For the v = 2 



comes better as N increases a conjecture already made by 
iMalyshkin & Goodman! (1200 lh . 

3.4. Fragmentation 

The grouping method can also be used in steady-state sit- 
uations where fragmentation balances coagulation. From the 
simulations in § l3.1l we expect that for the distribution method 
this is not a problem, as it follows both the number as well as 
the mass density. We will therefore concentrate our efforts 
on the equal-mass method. Fragmentation events in Monte 
Carlo may seem problematic, because catastrophic destruc- 
tion can result in the creation of many small particles that con- 
sume computational resources. However, due to the grouping 
method we now have a natural way to deal with fragmenta- 
tion. In this section we consider a test case of complete de- 
struction - meaning: the breakup of a particle into its smallest 
constituents (monomers) - that occurs once a mass threshold, 
'Wfrag, is exceeded. We then solve for the emergent, steady- 
state, mass distribution. Note that due to stochastic fluc- 
tuations the number distribution is never really steady-state. 
We will therefore average over many distributions during the 
same simulation run to obtain the 'steady-state' mass spec- 
trum. 

In Fig. [9] the results of the sum kernel, K,j - mi + nij, are 
given. In Fig.|9^ the fragmentation threshold is nif mg = 10 10 . 
The thick crosses denote the average over ~ 100 distribu- 
tions during the simulation after the initial distribution has 
relaxed to semi-steady state. The volume of the simula- 



simulations we used a number of N* = 60000, whereas in tion is stabilized when it incorporates N* = 10 species at 

*V — 2 x 10 12 » mfr a g. Figure|9p presents the result for an 
even larger value of mw = 10 20 . The same characteristic dis- 
tribution - resembling a bath tub - is found. These results sug- 
gest that, for the sum kernel with catastrophic fragmentation 
at a cut-off mass, most of the mass density becomes located 
near the end points. The system may then be approximated 



the v = 3 simulations we take N* = 300000. At v = 3, 
runaway occurs at a lower value of Mn and we are allowed 
to start with a higher number of groups. However, the dis- 
crepancy between the simulations at different values of N be- 
comes large in terms of CPU-time: the highest- N simulations 
take about half a day to finish. 

The results for the kernels K - {mimj) 2 and K - 
{mimjf are presented in Fig. [8] We extend the model of 
IMalyshkin & Goodman] d2001b to an unprecedented 10 160 
number of particles (about the square of the number of par- 
ticles in the universe). To speed up the simulations at high 
N we have decreased the mass-ratio parameter f e to 10" 3 
(which we checked to be still sufficiently accurate). The num- 
ber of simulations varies between 7V s i m = 100 at low N and 
A^sim = 10 at N = 10 160 . We stop the simulations when the ra- 
tio between the two largest masses becomes Ma /Mli ^ 10 4 , 
which is a measure of t R . By this time the simulation has ei- 
ther already finished or has reached it asymptotic limit. 

According to the ana lytic prediction (eq. 11201 . 
Malyshkin & Goodman 2001) the fs-log/V power-law is 
-1 and -2 for v = 2 and v = 3, respectively. In Fig. [8] the 
piecewise exponent is indicated (number above the lines), 
which shows a declining trend, except for the N - 10 160 
point. To check convergence, we have additionally performed 
five simulations at a higher numerical resolution for N - 10 80 
and N = 10 I6() (N* = 2 x 10 5 for the v = 2 simulations and 

N* = 2 x 10 6 for the v = 3 simulations). For N = 10 80 it 
was found that the results had converged, but at N — 10 160 
a reduced f« compared to the lower resolution simulations 
was found. This increased the local slope between N - 10 80 
and N - 10 160 (dotted lines, hardly distinguishable from 
the dashed lines, with the slope indicated below the line). 
Altogether, Fig. [8] suggests that the agreement between 
the numerical simulation and the monotrophic model be- 



by two states at m - mo ~ 1 and m - ntf ~ mf rag . The mass- 
loss of the monomer state due to collisions with the high-m 
state is mdm/dt\~ - HoHfWf and the gain due to fragmentation 
within the high mass bin becomes mdm/dt\ + = n 2 m 2 f . Balanc- 
ing these we get «o = fiftrif, i.e., equal mass densities for the 
monomer and high-m state. 

However, with other kernels the mass distribution due to 
fragmentation may also e merge as a pow er-law. To see this 
we adopt a simple model (Camacho 2001) in which the mass 
sweep-up, m, of a particle of mass m due to particles of lower 
mass m' < m is assumed to change only the mass in a contin- 
uous way. That is, 



= I 



dm'K(m, m')m'f(m'), 



(21) 



will shift particles of mass m by an amount Am - mAt. The 
total particle change per unit time in the interval [m, m + Am] 
becomes -A[mf(m)] due to gradients in m and the distribution 
function. The other side of the distribution then acts as a sink 
to the particles of mass m, i.e., R(m)f(m)Am is the number 
of particles in [m, m + Am] that is 'consumed' by the more 
massive particles, with R(m) 



R(m) 



*J 111 



dm'K(m, m)f(m) 7 



(22) 



the removal rate. For a steady state we must therefore have 
d(mf(m)) 



dm 



-f{m)R(m). 



(23) 
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Steady state distributions of the sum kernel. At m ■ 
10 20 and Nt 



10 



: "'frag particles are catastrophically disrupted into monomers of mass m = mo = 1. (A) mf™ = 10 
and N* = 10 4 ; (B) mf ra „ = 10 20 and N* = 2 X 10 4 . Distributions are averaged over various times during the steady-state process. Crosses denote this average 
over ~ 100 distributions. The spread is indicated by the grey shading. 



Assuming a power-law for /, / oc m a , equation ( 1231 can 
be solved for a by extending the limits of the integrals in 
equation (l2lT i and equation (l22l to and co, respectively 
(Camacho 2001). If |/?| < 1 the integrals converge and the 
resulting power-law exponent is a - -(3 +/3)/2. The physical 
reason for the scale-free solution is thus that, contrary to the 
linear kernel, the behavior of the distribution near the cut-off 
points is unimportant for the (local) distribution. Other works 
have confirmed this power-law exponent, (e.g., White 1982; 
Havakawa 1987). If [3 > 1 the integrals do not converge: 
i.e., collisions with the boundaries of the distribu tion d omi- 
nate and no power-law emerges. However, Klett d 19751) has 
shown that the/3 = 1 kernel Kjj = (m/m^) 1 / 2 , which therefore 
somewhat resembles the sum kernel, has a a - -2 power-law. 

To test whether a power-law emerges, we have performed 
fragmentation simulations for the constant kernel, K = 1, see 
Fig. [10] Thus, as /3 — we expect to end up with a slope 
of 1 /2 for the m 2 f(m) mass distribution. This trend is indeed 
observed in Fig. [TOh but only at large masses. The reason why 
an extension of the power-law to low masses is not seen, is 
the same as in Fig.|5^. Again, due to the equal-mass method, 
region of low mass density are not represented. Additionally, 
the injection rate at m = mo is episodic: at one instance many 
monomers of total mass ~ /«f rag enter the system, such that 
the mass density at this particular instance becomes very high 
(in comparison to the expected power-law distribution) for a 
very brief period of time. However, as coagulation is a non- 
linear process, i.e., dn^jdt &- —r&, these particles coagulate 
very quickly with each other. In other words, the episodic 
injection rate has the effect of speeding up the coagulation 
compared to a smooth rate and, when averaged over time, the 
true distribution is underestimated at places where the particle 
density strongly fluctuates with time. 

However, this deficiency of the equal-mass method does not 
affect the results at higher m: here the distribution is well re- 
solved. Additionally, the mean (time-averaged) influx rate of 
monomers, /, at m = mo is available, as well as the removal 
rate of the higher mass particles, R(m). A simple remedy to 
retrieve the low-m distribution becomes then available. In a 
separate simulation we have focused only on the low-m part 



of the mass spectrum until a cut-off mass m cut , and addition- 
ally include the fixed rates / and R, which were obtained from 
the previous simulation. Thus, R(m) is the probability that a 
particle of mass m collides with any particle of the 'frozen' 
large-m distribution that was previously obtained in Fig.fTOa. 
Finally, if particles manage to grow above m cut (i.e., avoid 
to collide with a high-m particle through R) they are also re- 
moved. In this way we effectively re-perform the same frag- 
mentation simulation, but in a reduced simulation volume, 
where the much smaller inter-collision time-increments At are 
appropriate for the low-m tail. In Fig.flOb this procedure is il- 
lustrated. Here, we put the cut at m =s 10 5 , recalculated the 
steady-state for small m, and 'glued' it to the high-m distribu- 
tion of m > m cut from Fig.fTOa. 

4. SUMMARY AND DISCUSSION 
4.1. Merits of the Method 

In the previous section we have applied the grouping 
method to a variety of situations and found satisfactory be- 
havior, both in terms of accuracy as in required computing 
power, despite the large dynamic range under consideration. 
The collisional evolutions of the analytic kernels were tested 
and agreed with the theoretical distributions. We also con- 
firmed the predictions made for a class of runaway kernels, 
which required us to model an unprecedented large number 
of particles. Finally, we applied the grouping method to sit- 
uations in which fragmentation operates and found excellent 
agreement with earlier findings (at least for the constant ker- 
nel; we are not aware of a prediction for the shape of the dis- 
tribution function in the sum-fragmentation kernel). 

Here we further discuss the (dis)advantages of both the 
equal mass and distribution methods and compare them with 
other studies. In short our conclusions are the following: the 
equal mass method is a true MC method, in the sense that no 
merging of the species array is required, while the mass den- 
sity of the distribution is resolved (following a characteristic 
mass). For non-runaway situations this is therefore the pre- 
ferred metho d. Other work along these lines can be found in 
Shima et alj d2007l) . These authors implement a superparticle 
method to model rain drop formation. These super-droplets 
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Fig. 10. — (A) Catastrophic fragmentation in the constant kernel, K = 1, with mf lag = 10 lu . At masses above m ~ 10 4 a power-law develops with a slope 
approximately equal to the predicted slope (a = -3/2; dashed line). At low-m the simulation lacks resolution and the power-law relation breaks down because of 
the episodic nature of the monomer injection. (B) This discrepancy is solved by performing a follow-up simulation for the low-mass part of the distribution only 
(m < m cut = 10 5 ) but including a removal rate due to collisions with large masses above m cut , R(m), and a monomer injection rate, /, that were obtained in (A). 
The distribution above m cu , is copied from (A). 

(like our groups) represent a multiple number of physical par- 
ticles with the same properties (like mass) and position. How- 



ever, their collisional setup is different. Where we allow many 
particl es of one group to collide with another, in Shima et al. 
(2007) there is at m ost one collision per particle. Also, in the 



IShima et all (120071) work the number of superdroplets is con- 
served and therefore (if we understand their picture correctly) 
the mean mass of the superdroplets. Thus, they always need 
a large number of superdroplets in order to resolve the high- 
m flu ctuations (see the discussion on Fig. 2c of lShimaetal.1 
2007). The biggest difference, therefore, is that our method 
allows a dynamical adjustme nt of the group sizes. Thus, we 
think the Shim aetall (120071) work is not capable of achiev- 
ing the high dynamic range of our grouping method. Still, we 
do note that their method achieves growth factors of at least 
10 6 in mass, and this may be sufficient for application to the 
atmospheric conditions of rain formation. 

With the distribution method we are able to trace well the 
fluctuations of the mass (or distribution) spectrum, by choos- 
ing the zoom factors such that there are a fixed number of 
groups for all exponentially spaced mass bins. This method 
is therefore well suited for runaway kernels. However, the 
consequence is a rapid proliferation of the N s parameter as 
many new species (which may be insignificant by, e.g., mass) 
are created. To limit N s , these species need to be merged, 
which somewhat violates the spirit of MC-methods as this 
includes averaging, not only over mass but also over struc- 
tural parameters. The merging feature renders th is method 
similar to WetheriU's discrete 'batches' technique (Wetherill 
1990). However, in our method structural parameters, like 
porosity and charge, can be included and there is consid- 
erable (useful) freedom in how the averaging can be per- 
formed. These structural parameters survive, in averaged 
form, the merging process in our distribution method, and 
can be fundamental in the process of coagulation. For ex- 
ample, fractal grain structure (porosity) and dust charging can 
lead to, respectively, strong particle grow th and even gelation 
dKonopka et alJl2005tlOrrnel et alJl2007l) . 



Both the Wetherill ( 1990) method and our approach give the 
correct behavior for the most massive and second most mas- 
sive p articles in the case of the product kernel (see Ina ba et all 
1999 for verification). Still, WetheriU's method has to deal 
with a number of free parameters: the spacing between the 
discrete batches and a parameter that determines the timestep, 
which must be fine-tuned. Our MC-method, although not free 
of them either, does have the advantage that the free param- 
eters correspond to physically-identifiable quantities, such as 
the mass-fraction parameter, f e . Another key advantage of 
our code is that the timestep, Af, follows naturally from the 
MC-code and is not required to be specified a priori. Espe- 
cially in the super-runaway kernels with v = 2, 3, growth can 
be very erratic: during the runaway-process there is a large 
amount of growth during a small time step. The success of 
WetheriU's approach lies in 'discretizing' the original (mean- 
field) Smoluchowski equation in an efficient way. Naturally, 
in Monte Carlo codes, discreteness is assured. 

The existence of two separate methods might appear to be a 
drawback. That is, runaway kernels (modeled with the dis- 
tribution method) and information on structural parameters 
(through the equal mass method) are difficult to treat simul- 
taneously, especially because one does not know a priori if 
a runaway situation will materialize. However, information 
about the kernel is 'hard- wired' into the Monte Carlo program 
and it is therefore clear at which point the situation is suscepti- 
ble to runaway growth. Then, these particles (the high-ra par- 
ticles) have to be further resolved, and this can be achieved 
by switching to the distribution method. Thus, we feel that 
it is possible to model runaway kernels with MC-methods, 
provided one implements a hybride switch that signals run- 
away growth through the effective mass scaling exponent that 
a kernel exhibits during the simulation. For example, consider 
the cross section of equation (O. Obviously, the mass scaling 
of cr, m v , becomes stronger than unity (v « 4/3) for suffi- 
ciently large masses. It does so in a smooth manner, allowing 
a hybride switch to be implemented. That is, as the program 
signals a large value of v for a particular kernel, it switches 
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from the equal mass to the distribution method, i.e., assigning 
lower zoom factors for the high-m particles and, if needed, 
merging structural parameters. This kernel ada ptation tech- 
nique is analogous to adaptive mesh refinement dWise & Abell 
2007), where the latter is grid based/spatial and we take a par- 
ticle/stochastic approach. Both numerical treatments require 
one to look carefully at the consequences of the switch, i.e., 
the impact of the merging prescription or the under-resolved 
part of the grid. 

The bottom-line on MC-methods thus becomes clear: there 
is a priori no universal way in which the correct, accurate re- 
sult can be guaranteed. Each situation requires interpretation 
and, if necessary, (re)adjustment of some aspects related to 
the grouping method; most notably through the choices for 
the grouping factors, i.e., the {z,}. However, MC-methods are 
flexible enough to do so in physically meaningful way. As 
such, they are a very powerful tool to describe any dynami- 
cal system in which particles interact in a manner that can be 
statistically represented. In fact, the techniques presented in 
this work allow one to treat a wide range of, astrophysically 
relevant, kernels. The nicest feature of MC-methods is that 
besides mass, structural properties of particles, like porosity 
and charge, can be directly incorporated into the collisional 
dynamics. 

4.2. Astrophysical implications 

There are several astrophysical situations where a large dy- 
namic range is required and for which our algorithm may 
be applicable. In protoplanetary disks, one of the key ques- 
tions is how the tiny (sub) micron-sized grains are transferred 
to planetesimals (Domin iket alJ 120071) . Because of the in- 
creasing relative velocities with particle growth, it may very 
well be the case that, while grains grow, the likelihood of 
fragmentation increases. An additional argument for frag- 
mentation is that collisional growth very efficiently removes 
the smallest grains, contrary to what o bservations indicate 
dDullemond & Dominikl2005l) . Recently, iBrauer et all d2007l) 
has solved the steady-state coagulation-fragmentation equa- 
tion for the dust distribution. It would be worthwhile to 
investigate whether these results are also obtained with the 
grouping method, and how this picture changes by inclu- 
sion of structural parameters, e.g., the porosity of aggr egates 
(IOssenkopfll993blKempf eTatlfl999HOrmel et alJl2007h . 

If particles reach the planetesimal size (~ 1 km) self-gravity 
may be sufficient to keep the collision products together. Fur- 



thermore, gravitational focussing, enhancing the gravitational 
cross-section, wil l turn the accretion pro cess into a runaway 
process (v > 1). Goldreich et al. (2004)) ha s outlined under 
which conditions runaway growth proceeds. iGreenberg et alJ 
( 1978) followed the collisional evolution (including fragmen- 
tation) of a population of initially 10 12 km-sized planetesimals 
until sizes of ~ 500 km were reached (after which particle-in- 
a-box approximations fail). We again note that a bin-based 
approach may have difficulty to model the runaway stage un- 
less it uses an ad-hoc extension to include the runaway particle 
explicitly (Wetherill 19901) . iBromlev & Kenvonl ([2006 ) have 
developed a hybrid code that combines a N-body and statisti- 
cal approach. However, these models can only treat one vari- 
able; with a Monte Carlo approach the internal structure of 
the particles can be included and a rubble pile structure can 
be distinguished from a homogeneous one. 

In stellar clusters, too, conditions may be feasible for run- 
away growth. Due to dynamical friction (equipartition) the 
more massive stars will s ink to the center of the cluster and 
we may even have v — 2 (Malyshkin & Goodman 2001; Lee 
2000). Moreover, the cluster will eject the smaller members 
and contract as a whole. These processes lead to core collapse 
in which the inner cut-off for the power law density function, 
ro, disappears, ro — » 0, and the density seemingly becomes 
singular. This condition is very favourable for a runaway- 
growth scenario, and may well be the mechanism for the for- 
mation of h igh mass stars, seed black holes or even gamma 
ray bursts dPortegies Zwart et alJ 119991; iFreitag et all 1 120061; 
iPortegies Zwart & van den Heu vel 2007). Globular clusters 
are already successfully modeled by MC-methods where the 
structural parameters are now the energy and angular momen- 
tum characterizing an orbit on which a particle (group) resides 
(e.g. iJoshi et aTll2000l) . Here, the grouping mechanism may 
then be well suited to model systems with a very large number 
of stars and stellar ejecta. Similarly, the collapse, and possi- 
ble fragmentation, of an interstellar gas cloud, leading to the 
formation of a stellar cluster or black hole, could be treated 
by following many (porous) gas clumps. 



We appreciate the efforts of Volker Ossenkopf to proof-read 
the manuscript and the useful suggestions and comments of 
the two referees, which significantly improved the structure 
and presentation of this work. 
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APPENDIX 
GROUPED COLLISION RATES 

We justify the expressions for the grouped collision rates, equation dlOl . where the collision rates are divided by the number 
of collisions during the group collisions. We show that, provided the individual rates C,j do not change as particles accrete, the 
mean collision rate is independent of the grouping. 

Consider a species of N identical particles that react with an individual rate «o with an external particle. In the language of the 
grouping method of §|2]the N particles are of species j that, during the group process, react with a particle from species i and 
«o = C/j. Let the process be such that the N particles react in groups of N/w, where 1 < w < N is the total number of such events. 

Let Pf be the probability that the system is in state z, meaning that the j-th event (out of w) has occurred. At t — the system is 
in the first state, so Po(t - 0) = 1 and P,(t = 0) = for i + 1, Let the rate at which the system evolves from state i to the next 
state (i + 1) be given by At, then 

at 

—± = -A 1 P l +A P ; (A2) 

at 

§ = -AiPt + Ai-iPi-x; (1 < » < w) (A3) 

at 

dP 

—Z = +A w . l P w . 1 . (A4) 

at 

Let R(t) be the mean reactivity of the system measuring the mean number of particles that has reacted at time t (0 < R < 1). 
Likewise, let R 2 measure the square number of reactants, i.e., 

If we weigh the state-equations (eq. 1A4ID with i/n and {i/n) 2 , respectively, and sum over i we obtain: 

dR 1 ^ 
dt w -f-J 

rIR 2 ~ 1 "' 

/=0 

Now the At are proportional to the amount of particles left, i.e., At oc (1 — i/w)N. In the case of group collisions, furthermore, 
we have argued that the rates should be divided by the number of reactants, such that the collision rate stays the same. We thus 
obtain 

Ai — (w — i)ao, (A7) 

where ao is the individual collision rate between two single particles. Equation ( IA6t then transform to 

dR i -\ 

— =a (l-R); (A8) 

dW l -^ 2w-l- 1 \ 

— -=a \-2R 2 + R + -\. (A9) 

dt \ ww] 

This first equation of the first moment (average value) is independent of w. Its solution is of course the well known exponential 
decay corresponding to Poisson statistics, i.e., 

R(i) = 1 - exp[-o/]. (A10) 

Inserting R into the equation for the second moment we solve for R 2 as 

¥(t) = (l - e- c ") (l - e- c " (l - uT 1 )) . (Al 1) 
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The spread <x in the distribution then becomes 

-at 

o- 2 (w) = (l -e- c ") . (A12) 



w 



That is, for w — » oo, cr — > and we are certain about the fraction of objects that have collided. For small w, however, there can be 
a significant spread. 

By choosing the collision rates as in equation (IA7t - which is simply the rate for a collision with one particle, ao(N - iN/w), 
divided by the number of particles in the group (N/w) - we end up with equation (lAlOt for the mean number of collisions. Since 
this is independent of w the mean reactivity of the system is also unaffected by the choice of w. For the mean reactivity it does 
not matter whether w - N (no grouping) or w — 1 (all particles in the same group) as long as ao stays constant during the 
group process. However, the spread in R does depend on the choice of w, which is quite natural, of course. A similar procedure 
cannot be performed for the in-group collisions since here equation ( !A4b involves r terms. The group collision rates for in-group 
collisions are not exact. However, the rates of equation (fT0l > are still a good approximation, and will become exact for large w. 
Besides, in-group collisions are a relatively minor occurrence. 



